

% This code is written by Dr. Hungtang Ko.
% For a manuscript titled "Collective movement of schooling fish reduces locomotor cost in turbulence".
% Last modified on 4/17/2024

close all

%% load semi-automatically tracked fish nose using DLTdv8 and parse the data
v = load('data/NoseTrackIndividual1@6BLsLaminar_dvProject.mat');
v = v.udExport.data.xypts;
yl = full(v(:,2));

v = load('data/NoseTrackIndividual1@6BLsTurbulent_dvProject.mat');
v = v.udExport.data.xypts;
yt = full(v(:,2));
yt(end) = [];

v = load('data/NoseTrackSchool2@6BLsLaminar_dvProject.mat');
v = v.udExport.data.xypts;
yl_school = [full(v(:,4))];

v = load('data/NoseTrackSchool1@6BLsTurbulent_dvProject.mat');
v = v.udExport.data.xypts;
yt_school = [full(v(:,2))];

% create a mask and turn the untracked points (defaulted to 0 in DLTdv8) to nans 
mask = (yt_school~=0); 
mask = double(mask);
mask(mask==0)=nan;
mask = mask(:,1);

% conversion from px to BL
yl = yl/277;
yt = yt/268;
yl_school = yl_school/177;
yt_school = yt_school/172;

% time
t = (0:1250)/125;

% plot time series
figure(1)
subplot(4,1,1)
plot(t,yl)
subplot(4,1,2)
plot(t,yt)
subplot(4,1,3)
plot(t,yl_school)
subplot(4,1,4)
plot(t,mask.*yt_school)

%% fft

% trim all the time series to 5 s to ensure fish is visiable for all conditions at all times 
yl = yl(t<5);
yt = yt(t<5);
yl_school = yl_school(t>2.5 & t<=7.5);
yt_school = yt_school(t>0.2 & t<=5.2);
t = t(t<5);

% set axis limit
xlim1 = [0,30];
ylim1 = [0,0.02];
xlim2 = [0,5];
ylim2 = [-0.5,0.9];

% fft
Fl = fft(yl);
Ft = fft(yt);
Fl_school = fft(yl_school,[],1);
Ft_school = fft(yt_school,[],1);

% frequency axis
fs = 125;
f = (0:length(Fl)-1)*fs/length(Fl);

% plotting
figure(2)
subplot(4,1,1)
stem(f,abs(Fl)/length(t)*2)
xlim(xlim1)
ylim(ylim1)

subplot(4,1,2)
stem(f,abs(Ft)/length(t)*2)
xlim(xlim1)
ylim(ylim1)

subplot(4,1,3)
stem(f,abs(Fl_school)/length(t)*2)
xlim(xlim1)
ylim(ylim1)

subplot(4,1,4)
stem(f,abs(Ft_school)/length(t)*2)
xlim(xlim1)
ylim(ylim1)

% generate excel with raw data
writematrix(["frequency(Hz)","laminar indivdual (BL)","turbulent individual (BL)","laminar school (BL)","turbulent school (BL)"],'nose_oscilation_freq.csv','WriteMode','overwrite')
writematrix([t',abs(Fl)/length(t)*2,abs(Ft)/length(t)*2,abs(Fl_school)/length(t)*2,abs(Ft_school)/length(t)*2],'nose_oscilation_freq.csv','WriteMode','append')

%% filtering low freuqency signals from high freuqency oscilations
cutoff = 3; % unit hz
Fl(f>cutoff & f<125-cutoff) = 0;
Ft(f>cutoff & f<125-cutoff) = 0;
Fl_school(f>cutoff & f<125-cutoff) = 0;
Ft_school(f>cutoff & f<125-cutoff) = 0;

figure(3)
subplot(4,1,1)
hold all
plot(t,yl-ifft(Fl))
xlim(xlim2)
ylim(ylim2)

subplot(4,1,2)
hold all
plot(t,yt-ifft(Ft))
xlim(xlim2)
ylim(ylim2)

subplot(4,1,3)
hold all
plot(t,(yl_school-ifft(Fl_school)))
xlim(xlim2)
ylim(ylim2)

subplot(4,1,4)
hold all
plot(t,(yt_school-ifft(Ft_school)))
xlim(xlim2)
ylim(ylim2)

% generate excel with filtered data
writematrix(["t(s)","laminar indivdual (BL)","turbulent individual (BL)","laminar school (BL)","turbulent school (BL)"],'nose_oscilation_filtered.csv','WriteMode','overwrite')
writematrix([t',yl-ifft(Fl), ...
            yt-ifft(Ft), ...
            yl_school-ifft(Fl_school), ...
            yt_school-ifft(Ft_school)],'nose_oscilation_filtered.csv','WriteMode','append')



